Meso-scale computer modeling of lipid-DNA complexes for gene therapy 
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We report on a molecular simulation method which captures the self-assembly of cationic lipid- 
DNA (CL-DNA) gene delivery complexes. Computational efficiency required for large length- and 
time-scale simulations is achieved through a coarse-grained representation of the intra- molecular 
details, and via inter- molecular potentials, which effectively mimic the hydrophobic effect without 
explicit solvent. In addition to showing spontaneous self-assembly of complexes, the broad utility of 
the model is illustrated by demonstrating excellent agreement with X-ray diffraction experimental 
data for the dependence of the spacing between DNA chains on the concentration of CLs. At high 
concentrations, the large electrostatic pressure induce the formation of pores in the membranes 
through which the DNA molecules may escape the complex. We relate this observation to the origin 
of recently observed enhanced transfection efficiency of lamellar CL-DNA complexes at high charge 
densities. 
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Modeling and simulation of biological systems possess 
significant challenges due to the spatial complexity of 
such systems and by the accompanying range of tem- 
poral scales involved. In the case of lipid bilayers, for 
instance, the width of the bilayer is in the nanometer 
range, whereas the lateral dimension may extend up to 
several micrometers. The corresponding time scales span 
an even larger range of orders of magnitude. To address 
the multi-scale nature of membranes, a variety of models 
have been devised which differ in the length- and time- 
scale of the phenomena of interest pj. The two most 
extreme approaches are: (1) atomistic computer simula- 
tions in which the lipids and the embedding solvent are 
modeled explicitly in full detail 0], and (2) phenomeno- 
logical models such as the effective Helfrich Hamiltonian 
where the bilayer membrane is treated as a smooth con- 
tinuous surface with a small number of elastic coefficients 

□ 

The gap between those two approaches can be bridged 
by coarse-grained (CG) computer models in which groups 
of several atoms are represented by unified particles in- 
teracting via effective potentials 4] . The essential benefit 
of such models is that they have condensed the most de- 
tailed variations in both time and space, while retaining 
the fundamental properties of membrane self-assembly 
and elasticity. The level of detail employed by different 
CG models varies greatly. Some CG models are con- 
structed to mimic specific lipidic systems by selecting 
simplified representations for water molecules and dif- 
ferent chemical groups which constitute the lipids, and 
by developing effective force-fields that reproduce key 
structural features known from experiments and atom- 
istic simulations p. In coarser-grained models the am- 
phiphilic molecules are constructed from only two types 
of particles - hydrophilic and hydrophobic - which, re- 
spectively, are attracted and repelled from a third kind 



of particles representing the solvent 6] . Such models ad- 
dress more general features of self-assembly rather than 
specific systems. 

More recently, a number of research groups have devel- 
oped CG models in which the bilayer membranes are sim- 
ulated without direct representation of an embedding sol- 
vent; this is accomplished by constructing intermolecular 
force fields that mimic effects of hydration 0, H, ©, EU . 
The development of these "water-free" models is an im- 
portant accomplishment in large-scale membrane simu- 
lations, considering the fact that the number of solvent 
particles in conventional CG models is significantly larger 
than the number of lipids. Most of these models em ploy 
simple Lennard- Jones (LJ) type pair-potentials [3.l^.ll0j . 
and exhibit spontaneous self-assembly to a membrane as 
demonstrated in Fig. [I] (A). 

The computational simplicity and efficiency of water- 
free CG models make them excellent tools for testing 
mesoscopic scale theories of membranes and gaining in- 
sight into the basic molecular mechanisms that govern 
these systems. More importantly, these models may serve 
as platforms for CG simulations of more complex systems 
including bio-polymers and networks. This paper reports 
on such a computational study of complexes of cationic 
lipid and DNA (CL-DNA) in the context of recent exper- 
iments on molecular assemblies. The molecular assem- 
blies are formed spontaneously when DNA is mixed with 
cationic and neutral lipids in an aqueous environment 
[Til . Their formation is driven by the electrostatic 
attraction between the negatively charged DNA and the 
cationic lipid head-groups, and by the entropic gain as- 
sociated with the concurrent release of the ti ght l y b ound 
counterions from the CL and the DNA 111 13- x ~ 
ray diffraction experiments have revealed that CL-DNA 
complexes exist in a variety of mesoscopic structures > 
including a lamellar phase where the DNA is intercalated 
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FIG. 1: (A) Self-assembly of a bilayer patch of 100 lipids. 
Lipids are modeled as trimers where the black and grey 
spheres represent hydrophilic and hydrophobic particles, re- 
spectively (see details in Refs. Each Monte Carlo time 
step includes (on average) an attempt to move and rotate 
each molecule. The membrane in the last figure has been ro- 
tated for a better viewing of its structure. (B) Self-assembly 
of a CL-DNA complex consisting of a total number of 240 
lipids (of which 90 are charged). Red rods represent DNA 
segments. (C) Equilibrium configuration of a complex with 
cj) c ~ 0.9 whose membranes develop pores. 

in between lipid bilayers 11]. The DNA strands form 
a one-dimensional (ID) ordered array, where the DNA 
interaxial spacing c^dna decreases with the charge den- 
sity of the membranes. Isoelectric complexes, where the 
charges on the DNA exactly match those on the CL, are 
the most stable ones since they enable nearly complete 
counterion release . 

CL-DNA complexes have attracted much attention be- 
cause of their potential use as non- viral transfection vec- 
tors in gene therapy 0, EH • Cationic liposome 
transfer vectors exhibit low toxicity, nonimmunogenicity, 
and ease of production [TEL H^ . but their transfection 
efficiency (TE) remains low compared to that of viral 
vectors [l^ |2(| . This has spurred an intense research ac- 
tivity aimed at enhancing TE 0, 0, [2(J . Recognizing 
that the structure of CL-DNA complexes may strongly 
influence their function and TE, much of the effort in 
theoretical and experimental studies has been devoted to 
understanding the mechanisms governin g co mplex forma- 
tion, structure, and phase behavior fl3 l PuL l2ll l22l l23q . 

To the best of our knowledge, only a single attempt 
has so far been made to study CL-DNA complexes us- 
ing molecular simulations |24|. That study employed 
a fully atomistic description of both the lipids and the 



DNA and has, consequently, been limited to a small sys- 
tem consisting of 48 lipids (of which 20 were charged) 
and a short DNA segment of 10 base pairs. The dura- 
tion of the simulations was of several nanoseconds. The 
mesoscopic regime encompassing the statistics and evo- 
lution of large molecular ensembles is, however, inacces- 
sible to full atomistic simulations due to the computa- 
tional requirements, and only continuum behavior of ex- 
isting CL-DNA structures can be addressed based on free 
energy functionals which are insensitive to the fine de- 
tails of the lipids and DNA [2]|. Thus, our understand- 
ing of many important CL-DNA complex features, such 
as the process of self-assembly and membrane evolution, 
mesoscopic structure and defects, can be addressed only 
through CG simulations. 

The model presented here is based on the CG water- 
free membrane model, which produces self-assembled bi- 
layer membranes such as the one depicted in Fig. D](A); 
model details are given in Refs. |8|- The lipids consist 
of one hydrophilic (representing the head group) and 
two hydrophobic (representing the tail) spherical par- 
ticles with short-range pair-interactions between them. 
We set the diameter of these particles a 6. 3 A (see 
the definition of a in Ref. 8]). For this value of a the 
area per lipid in the original model au P id — 70 A 2 . In 
the present model a fraction </> c of the hydrophilic head 
groups carry charge +e. DNA is modeled as a rigid rod 
with a uniform axial charge density Adna = — e/1.7A and 
radius Rbna = 10A. Excluded volume interactions be- 
tween rods (R) and spheres (S) are introduced via a trun- 
cated and shifted potential of the form: UrsI^bT = 

50{[(a/2 + iW)/r] 12 -l}, where r is the distance 

between the center of the sphere and the axis of symme- 
try of the rod. The distance between nearest-neighbor 
rods is restricted to g^dna > 2_Rdna- 

We study isoelectric complexes where the total charges 
of the DNA and the CLs neutralize each other, with 
no added counterions. Simulations of the quasi two- 
dimensional (2D) complex are conducted in a rectangular 
system of size L x x L y x L z , with full periodic bound- 
aries along the x and y directions [25j |. and periodicity 
with respect to only lipid mobility and short-range inter- 
actions in the z direction. Simulations were performed 
at room temperature and with a bulk water uniform di- 
electric constant e r = 78. The rods are arranged in a 
one-dimensional array with equal spacing in the x?/-plane 
along the ^/-direction. We first examined the ability of 
the complex to self-assemble spontaneously: 240 lipids, 
of which 90 are charged, were randomly distributed in 
a box of size L x = 94A, L y = 5lA, L z = 157A. Three 
equally-spaced DNA strands, each carrying a total charge 
— 30e, were placed at the mid-plane z = L z /2 of the box. 
Canonical ensemble Monte Carlo (MC) simulations were 
used to generate the temporal evolution of the lipids, 
while the positions of the rods were fixed. Each MC 
step consists of (on average) an attempt to translate (and 
make some minute changes in the relative locations of the 
three particles with respect to each other) and rotate each 
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lipid. Typical evolution of the system is depicted in Fig. [I] 
(B). The Figure demonstrates that a complex is formed 
quite rapidly, but full association of the lipids is much 
slower and has not been accomplished by the end of the 
MC run. The large majority of CLs (except for, typically, 
as few as 3-4 molecules) tend to reside in the two inner 
monolayers facing the DNA array, as expected in view of 
the strong electrostatic attraction between the CLs and 
DNA. It is possible that the final snapshot in Fig. ^ (B) 
is typical for an equilibrium distribution of free and com- 
plexed lipids, although long MC runs of pre-assembled 
complexes showed no evidence for escape of lipids from 
the complex. 

One of the parameters which is believed to strongly 
influence the TE is the membrane charge density gm = 
e0 c /anpid |2fl[. It was found in X-ray diffraction experi- 
ments that for isoelectric complexes the DNA interaxial 
spacing, g^dna, is related to gm by [Til 
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Measuring g^dna and verifying this result can serve as a 
test for our model's ability to mimic the meso-scale be- 
havior of CL-DNA complexes. We have therefore per- 
formed a set of MC simulations of pre-assembled iso- 
electric complexes consisting of 300 charged lipids, an 
array of five equally spaced DNA molecules of length 
L y = 102 A (i.e., each carrying a total charge — 60e), and 
an additional number of neutral lipids. Simulations have 
been performed using the constant surface tension ensem- 
ble (iV, 7, T), with 7 = 09, 2^. The area of the complex 
has been changed by allowing flexibility in L x (rescaling 
the x coordinates of particles accordingly), keeping L y 
fixed. Attempts to change the area of the system were 
made, on average, twice every MC time unit. Two flat 
bilayers, each consisting of 150 cationic and (2N — 150) 
neutral lipids, were initially set up on both side of the 
DNA array with the head groups in the inner (outer) 
monolayers at distance 0.5g -\- Rbna (5.5<j + ^dna) from 
the complex mid-plane. Each monolayer consisted of N 
lipids, with the CLs occupying the two inner monolayers 
while the outer monolayers consisting of neutral lipids 
only. This distribution of lipids remained unchanged over 
the course of the simulations (i.e., there was no diffusion 
of lipids between the monolayers on the time scale of 
the simulations). The fraction of charged lipids is de- 
fined as the number ratio of charged to neutral lipids in 
the inner monolayers, <p c = 150/iV. The pre-assembled 
complexes were allowed to equilibrate for 5 x 10 5 steps 
prior to data collections over, at least, 2 x 10 6 additional 
steps. Membrane fluidity was verified by lateral diffusion 
of both neutral and charged lipids. The former diffuse 
faster than the latter. 

The average spacing between adjacent DNA rods, 
g^dna i is plotted in Fig. as a function of the inverse 
of the fraction of charged lipids, 1/0 C . The numerical 
data is in excellent agreement with the experimental re- 
sults reported in Ref. pj|. The solid line is obtained 



70 



60 



50 



o Q 40 



30 



20 



1% 



FIG. 2: Average DNA spacing, c?dna as a function of the 
inverse of the fraction of charged lipids l/<j> c . Circles - nu- 
merical results (the uncertainties in the data are smaller than 
the symbols); solid line - fit to Eq. JIJ with an p id = 69 A 2 . 



from Eq. JIJ with an p id = 69 A 2 . The deviation from 
linear behavior at high charge densities arises from the 
increase in the area per lipid as (j) c (see Fig. [3 left 
?/-axis, black circle symbols). Similar deviation is ob- 
served in the experiment |27j . The assumption under- 
lying Eq. is that the effective interactions between 
the DNA are repulsive and balanced by the elastic mem- 
brane forces. According to Hooke's law, the elastic stress 
acting on a membrane is related to an p id and its equi- 
librium value a° ipid by r = K A (a lipid - a° ipid ) /ag pid , 
where Ka is the 2D stretching modulus, which for lipid 
bilayers is typically in the range Ka ^ 10 2 ergs/cm 2 . 
At high charge densities, the electrostatic stress is suf- 
ficiently large to eliminate the membrane thermal un- 
dulations and increase aii p id [2q |. In the present study, 
we find for complexes with <p c ~ 0.85 that the strain 

£ = (aiipid - a£ p id) / a n P id ~ O- 1 ' whicn is the typical 
strain lipid membranes can withstand before rapture |2<j . 
Membranes with higher <j) c have indeed been found to 
be susceptible to pore formation, as illustrated by the 
configuration in Fig. ^ (C) of a complex with <j) c ~ 0.9. 
The loss of mechanical stability is also evident from the 
rapid decrease in the complex stretching modulus K\ for 
4> c ^ 0.7 (Fig. 01 right y-axis, red circle symbols), which 
has been extracted from the mean square fluctuation in 
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larger area fluctuations at high <p c increase the probabil- 
ity of pore opening which in turn may lead to disassoci- 
ation of the complex. We suggest that this may be the 
origin of the recently observed enhanced TE of lamellar 
CL-DNA complexes at high charge densities |20| . Trans- 
fection is viewed as a two-stage process: (1) cellular up- 
take via endocytosis, and (2) escape of the complex from 
the endosome, presumably through fusion of the lipids 
with the endosomal membrane and release of the DNA 
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FIG. 3: (a) Left y-axis, black circle symbols - average area 
per lipid an p id as a function of the fraction of charged lipids 
(j) c - (b) Right 2/-axis, red square symbols - Stretching modulus 
K\ of the complex as a function of <j) c . 



into the cytoplasm. TE of lamellar complexes is limited 
by the rate of the second stage and, hence, increases with 
the decrease of mechanical stability, i.e., with increase of 
charge density. 

In summary, the structural-mechanical properties of 
CL-DNA complexes have been studied using a coarse- 
grained "water- free" model. Despite of the highly sim- 
plified representation of lipids and DNA, the model re- 



liably demonstrates self-assembly of both a charge neu- 
tral bilayer lipid membrane and a lipid membrane DNA 
complex. The model reproduces accurately the mea- 
sured interaxial spacing c^dna over a wide range of charge 
densities. The wide applicability of the model should 
be attributed to the fact that the meso-scale behavior 
of the system is dominated by non-specific electrostatic 
and elastic interactions. At high charge densities we ob- 
serve that the increasing electrostatic pressure exerted 
on the membranes leads to pore opening through which 
the DNA may be released from the complex. The DNA 
release to the cytoplasm is a necessary step for successful 
transfection and, thus, our results explain the enhance- 
ment of transfection efficiency at high concentrations of 
CLs. Given the consistency of agreement between our 
course-grained molecular approach and observed exper- 
imental features, we suggest that the presented model 
is an appropriate and promising tool for investigating 
the statistics and dynamics of lipid-DNA complexes on 
spatial and temporal scales relevant for biological and 
biomedical applications. 

We are grateful to C. R. Safinya for many useful dis- 
cussions. OF and PP acknowledge the support of the 
MRL Program of the National Science Foundation under 
Award No. DMR00- 80034 and NSF Grant No. DMR02- 
037555. 



[1] M. Miiller, K. Katsov, and M. Schick, J. Polymer Sci. B: 
Polymer Phys. 41, 1441 (2003), and Refs. therein. 

[2] L. Saiz, S. Bandyopadhyay, and M. L. Klein, Bioscience 
Rep. 22, 151 (2002), and Refs. therein. 

[3] R. Lipowsky and E. Sackmann Eds. Structure and Dy- 
namics of Membranes (Elsevier, Amsterdam, 1995). 

[4] S. O. Nielsen, C. F. Lopez, G. Srinivas, and M. L. Klein, 
J. Phys.: Condens. Matter 16, R481 (2004). 

[5] E. Lindahl and O. Edholn, J. Chem. Phys. 113, 3882 

(2000) ; S. J. Marrink and A. E. Mark, J. Phys. Chem. B 
105, 6122 (2001). 

[6] R. Goetz and R. Lipowsky, J. Chem. Phys. 108, 7397 
(1998). 

[7] H. Noguchi and M. Takasu, Phys. Rev. E 64, 041913 

(2001) . 

[8] O. Farago, J. Chem. Phys. 119, 596 (2003); O. Farago 

and P. Pincus, J. Chem. Phys. 120, 2934 (2004). 
[9] G. Brannigan, A. C. Tamboli, and F. L. H. Brown, J. 

Chem. Phys. 121, 3259 (2004). 
[10] L R. Cooke, K. Kremer, and M. Deserno, preprint 

|cond-m at/0502418l 
[11] J. O. Radler, I. Koltover, T. Salditt, and C. R. Safinya, 

Science 275, 810 (1997). 
[12] J. O. Radler, I. Koltover, T. Salditt, and C. R. Safinya, 

Science 281, 78 (1998). 
[13] D. Harries, S. May, W. M. Gelbart, and A. Ben-Shaul, 

Biophys. J. 75, 159 (1998). 
[14] C. R. Safinya, Curr. Opin. Struct. Biol. 11, 440 (2001). 
[15] P. L. Feigner, Sci. Am. 276, 86 (1997). 
[16] P R. Clark and E. M. Hersh, Curr. Opin. Mol. Ther. 1, 



158 (1999). 

[17] S. Li and L. Huang, Gene Ther. 7, 31 (2000). 

[18] H. Kamiya, H. Tsuchiya, J. Yamazaki, and H. Ha- 

rashima, Adv. Drug Delivery Rev. 52, 153 (2001). 
[19] H. F. Willard, Science 290, 1308 (2000). 
[20] K. Evert et a/., Curr. Med. Chem. 11, 133 (2004). 
[21] S. May and A. Ben-Shaul, Curr. Med. Chem. 11, 151 

(2004). 

[22] R. Bruinsma, Euro. Phys. J. B4, 75 (1998). 

[23] D. Harries, S. May, and A. Ben-Shaul, J. Phys. Chem. B 

107, 3624 (2003). 
[24] S. Bandyopadhyay, M. Tarek, and M. L. Klein, J. Phys. 

Chem. B 103, 10075 (1999). 
[25] N. Gr0nbech- Jensen, G. Hummer, and K. M. Beardmore, 

Mol. Phys. 92, 941 (1997). 
[26] An alternative scheme [O. Farago, in preparation (2005)] 

for sampling the (iV, 7, T) ensemble has been used in this 

work. 

[27] The experimental results reported in [see Fig. 4 (B)] 
show agreement with Eq. (J at low charge densities 
(with a value of an p id which is slightly different than 
the one defined here) and a similar deviation trend at 
high charge densities. Note that our Eq.Jfl) and the com- 
parable one in Ref . [TT| express the same relationship in 
different forms. See discussion in I. Koltover, T. Salditt, 
and C. R. Safinya, Biophys. J. 77, 915 (1999). 

[28] A. W. C. Lau, and P. Pinucs, Eur. Phys. J. B 10, 175 
(1999). 

[29] D. Boal, Mechanics of the Cell (Cambridge Univ. Press, 
Cambridge, 2002). 



